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Abstract. GTPase molecules are important regulators in cells that continuously run through 

an activation/deactivation and membrane-attachment/membrane-detachment cycle. Activated 

GTPase is able to localize in parts of the membranes and to induce cell polarity. As feedback loops 

I contribute to the GTPase cycle and as the coupling between membrane-bound and cytoplasmic 

I processes introduces different diffusion coefficients a Turing mechanism is a natural candidate for 

f^ this symmetry breaking. We formulate a mathematical model that couples a reaction-diffusion 

^SJ system in the inner volume to a reaction-diffusion system on the membrane via a flux condition 

. .J and an attachment/detachment law at the membrane. We present a reduction to a simpler non- 

flj local reaction-diffusion model and perform a stability analysis and numerical simulations for this 

/**N reduction. Our model in principle does support Turing instabilities but only if the lateral diffusion 

^^ of inactivated GTPase is much faster than the diffusion of activated GTPase. 

1^ 

r-— I 1. Introduction 

-^ GTP-binding proteins (GTPases) are crucially involved in many processes in cells such as 

membrane traffic, cellular transport, signal transduction, or cytoskeleton organization [301 I12j . 
Common to the diverse families of GTPase is the cycling between an active and an inactive 
state. Besides the activation-inactivation cycle there is also a spatial cycle: in the cytosol almost 
all GTPase is inactive whereas the active state is only present at the membrane. Reaction 
and diffusion processes both in the cytosolic volume and on the membrane surfaces as well as 

^N membrane attachment and detachment contribute to the proper function of GTPase molecules. 

,^ For different GTPase localization into subcellular compartments has been observed and has 

Q^ been recognized as crucial for its function. Cluster formation of activated small GTPase Cdc42 

•/"J precedes the budding of yeast p'5j , other small GTPase of the Rho-subfamily are known to form 

micro-domains on continuous membranes \'27\ [^ . Such a transition from a homogeneous distri- 

^^ bution to a polarized state is often key for the formation and maintenance of complex structures. 

^Z, The emergence of localized structures is typically driven by a continuous input of energy 

> 



B 



Turing |3H [23] pioneered models for symmetry breaking by diffusion-driven instabilities. These 
are based on a slowly diffusing self-activator and a highly diffusive antagonist |18j . Self-activation 
is typically present by some kind of feedback. In activator-substrate-depletion type Turing mech- 
rS anisms the production of the activator induces a decrease of the substrate. Diffusion-driven 

c3 instabilities typically require large differences in the diffusion coefficients of the activator and its 

antagonist. In many biological applications this is not realistic and Turing type mechanisms can 
therefore not explain symmetry breaking events. In our context, however, cytosolic diffusion is 
typically much faster than lateral diffusion. Coupled systems of 2D and 3D reaction-diffusion 
processes might therefore be a candidate for a realistic Turing mechanism. 

Distinct mathematical models for the GTPase cycle have been proposed and analyzed, with 
diverse conclusions. A general model for signaling molecules in a cell that cycle between a non- 
recruiting cytosolic state and a recruiting membrane-bound one has been evaluated in [1]. There 
the emergence of cell polarity has been demonstrated for an intrinsically stochastic mechanism 
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for self-activation by positive feedback. A corresponding deterministic model in contrast has 
been shown not to produce any heterogeneous pattern. However, the deterministic PDE model 
in [1] does not directly treat processes in the cytosol as all variables are membrane bound and 
all reactions are local. A complex PDE model that accounts for chemical reactions, membrane- 
cytoplasm exchange and diffusion is given in [9]. There a scaling factor accounts for the different 
volume of a thin 3D membrane layer and the inner volume. Variables representing membrane 
bound molecules and variables representing cytosolic quantities however both have the same 
domain of definition. Numerical simulations and a linear stability analysis show that the model 
allows for a Turing mechanism. The formation of micro-domains in a GTPase cycle model are 
also studied in [3j. Here the equations are formulated on a flat membrane surface. It is shown 
that no Turing pattern can occur unless an extra flux term is included. This flux term accounts 
for interactions between GTPase and membrane proteins and represents a phase separation type 
energy gradient. As an alternative explanation for the emergence of cell polarity in GTPase 
mediated processes a 'wave-pinning' mechanism is proposed in [22]. A two component system 
for the nucleotide cycle is suggested with a Hills type non-linearity that leads to a bistability. 
Domains are formed by emerging traveling waves that are stopped by a decreased supply of 
non-activated GTPase. 

Our goal is to introduce a model for the GTPase cycle with an improved coupling of processes 
with different dimensionalities. We investigate whether a Turing type instability - of activator- 
substrate depletion type - could possibly explain the localization of activated GTPase on the 
membrane. In Section [2] we will first formulate our mathematical model and derive a reduction 
that only incorporates membrane-bound active and membrane-bound inactive GTPase. We per- 
form a stability analysis and numerical simulations for this reduction. The explicit dimensional 
coupling in the full model is still reflected by the appearance of a non-local term. We show in 
Section [3] that for this model Turing patterns are possible. Our numerical simulations in Section 
[5] confirm this result and shed some light on the kind of patterns that are supported by our model 
and the influence of changes in different parameters. We develop here a general numerical scheme 
that can be extended to general membrane geometries and to more involved coupling laws. In 
Section |4] we investigate - even for a more general class of similar models - whether for equal 
diffusion constants of activator and substrate, i.e. activated and non-activated GTPase, Turing 
pattern are possible. Our results will finally be discussed in Section [6] 
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2. Model Description 

2.1. Mechanistic description of the GTPase cycle. Here we briefly review the key steps of 
the GTPase cycle as indicated in Fig. [T| Chemically, the difference between the active and inactive 
state of the GTPase is that in the active state guanine-tri-phosphate (GTP) is bound whereas 
guanine- di-phosphate (GDP) is bound in the inactive state. Only activated GTPase interacts with 
downstream effectors. Activation of a GTPase is by exchange of GDP by GTP, inactivation by 
hydrolysis and dephosphorylation of GTP to GDP. Both processes are intrinsically very slow and 
need the catalyzation by a GEF [guanine exchange factor) and GAP protein [GTPase activating 
protein), respectively |11|,[2]. Cytosolic GTPase can only be found in complex with a displacement 
inhibitor (GDI) that prevents the binding of GTPase to the membrane. As the affinity of GTPase 
towards GDI is much higher when GDP is bound, predominantly the inactive state occurs in the 
cytosol [8j . How GDP-bound GTPase is released from the complex with GDI and how it associates 
to the membrane is less clear, mediation by a GDI displacement factor (GDF) has been proposed 
as a possible mechanism [26J. For several GTPase positive feedback loops have been identified 
that support the activation of GTPase. Activated GTP-Rab5 is known to recruit a cytosolic 
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GEF-effector complex (RabexS and RabaptinS) to the membrane and to increase the activity of 
the GEF [IQj. A similar feedback loop has been found for activated Cdc42 GTPase [33]. In the 
following we formulate a mathematical model that reflects the key features of a GTPase cycle. 
Our main focus is on the treatment of the dimensional coupling and less on a detailed description 
of the reaction kinetics. 



GDI GDP , 





GEF 



Figure 1. The GTPase reaction cycle: The activation of GDP-bound GT- 
Pase is either catalyzed by GEF (lower semi circle) or by an effector-GEF-GTP- 
GTPase complex (upper semi circle). The inactivation of GTP-bound GTPase 
is catalyzed by GAP. Further reactions that are depicted are GDP-GDPase-GDI 
complex formation/dissociation and effector-GEF-GTP-GTPase complex forma- 
tion/dissociation. See the text for additional information. 



2.2. The mathematical model. We restrict ourselves to the investigations of processes in the 
cytosol and at the outer plasma membrane only. Inner organelles with additional membrane 
boundaries could be included as well. The cytosolic volume of a cell is represented by a bounded, 
connected, open domain i? C M^ and the cell membrane by the boundary of B that we assume 
to be given by a smooth, closed two-dimensional surface P := dB without boundary. In addition 
we fix a time interval of observation / := [0, T] C M. We formulate a system of PDE's for the 
following unknowns: 

y : i? X / — )• R concentration of cytosolic GDP-GTPase (in complex with GDI), 

V : r X / — ;■ M concentration of membrane-bound GDP-GTPase, 

u : r X / — ;■ M concentration of membrane-bound GTP-GTPase, 

m : r X / — )• R concentration of membrane-bound, effector-GEF-GTP-GTPase complex, 

^r : r X / — )■ R concentration of membrane-bound GEF. 
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We prescribe initial conditions at time t = 0, 

V{-,0) = Vo, v{-,0) = vo, u(-,0) = no, m(-,0) = mo, g{-,0) = go, 
Vq : B ^ R, vo,uo,mo,go : T -^ M. 

The coupling condition between cytosolic and membrane processes involves a Neumann boundary 
condition for V that is specified below. Physical units are given by 

IV] = '^. |„l = |„l = H = M = =^. 

m'^ m^ 

Much more extended sets of variables could be considered here. In particular we do not explicitly 
take into account the effector and GAP concentrations. Catalyzation of the inactivation process 
will be described implicitly. 

2.2.1. Reaction Kinetics. We assume simple mass action kinetics or a Michaelis-Menten type law 
for catalyzed reactions. For the change of concentration of the above variables due to reactions 
we prescribe the following equations. Our choices here are similar to the more general model in 

The concentration of membrane-bound GDP-GTPase is decreased by the activation process, 
which is catalyzed by both GEF and the effector-GEF-GTP-GTPase complex. For the corre- 
sponding rates we assume that they are proportional to the GDP-GTPase concentration and 
the concentrations of the catalysts. Vice versa GDP-GTPase is produced by the inactivation of 
GTP-GTPase. Since we have not taken the GAP concentration into account, we here assume a 
Michaelis-Menten law for the kinetics. The change of v due to activation and inactivation we 
therefore describe by 

u 
[9tw]rcaction = -kivg - k2vm + ks — — - on r X /. (2.1) 

U + K4 

For the change in u we have in addition to the processes above the production of u by dissociation 
of the effector - GEF - GTP-GTPase complex and the loss due to the formation of this complex. 
We model this by the equation 

[<9in] reaction = hvg + k2vm - k^ — — k^ug + k-^m on F x /. (2.2) 

u + ki 

Correspondingly complex formation and dissociation lead to the following laws for the concen- 
tration of the complex and GEF (or rather a GEF-effector complex as we do not take explicitly 
into account the effector). 



[9fm]reaction = k^ug - k^^m on r X /, 






(2.3) 


[5*5] reaction = -k5ug + k^^m OU F X /. 






(2.4) 


Units for the reaction rates ki are given by 
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[ki\ - M - M - ^^1 . ^ , M - ^2s ' ik,\ - ^2 , 





2.2.2. Simplified Kinetics. For later use in the mathematical analysis we use a quasi-steady state 
approximation for the complex formation to do a first reduction of our kinetic laws. We assume 

m = - — ug (2-5) 

fc_5 

and use GEF-conservation 

m + g = const =: ^o, (2-6) 



TURING INSTABILITIES 
mol 



where [go] = [g] = ^^. If we take initial data mo = we have go = go. Equations (2.5), (2.6) 
yield 

K^ugo 



m 



5 = 50 1 



l + K^u' 

K5U 



I + K5U 



with K^ := Y^ ciiid [K^] = ^^. From this, one obtains simplified rate equations for the change 
due to reactions 

\a } I, ( ^ ^5^ ^ 1 K^ugo u 

l9tvUcUon = -hvgo[l-^^^)-hv^^^ + ks^^^^^ on T x /, 

[t^ftt] reaction = hvgo ( 1 " -, , L ) + ^2^ -, ,^ r^" ^3 "l" = -[^tu] reaction On T X /. 

2.2.3. Diffusion. We describe cytosolic diffusion of the (inactive) GTPase in B by the standard 
Laplace diffusion operator A and a diffusion constant D > 0. Lateral diffusion on the membrane 
is described by the Laplace-Beltrami-operator Ap (which is the generalization of the ordinary 
Laplacian to surfaces [5] ) and diffusion constants du and dy for the active and inactive membrane- 
bound GTPase concentrations, respectively. 

[dtVUmnsion = DAV in B x I, (2.7) 

[5tM] diffusion = duAru on F x /, (2.8) 

[9tu]diffusion = d^Arv on F x /. (2.9) 

2.2.4. Membrane attachment and detachment. We describe the association and dissociation of 
inactive GTPase at the membrane by a flux boundary condition for F at F 

-DVV-u = q on F, (2.10) 

where v denotes the outer normal to B at F. For the flux q we formulate a constitutive equation: 
membrane attachment is treated as a reaction between cytosolic GTPase and a free site on the 
membrane and modeled by a Langmuir rate law [T7|. Detachment is taken proportional to the 
inactive GTPase concentration, which together gives the equation 

I Rl 
q = hQ---V{cra3.^-u-v)+-b-QV. (2.11) 

Here Cmax denotes a saturation value and 66,6-6 are sorption coefficients. By (cmax — u — v)+ 
we denote the positive part of Cmax — u — v as adsorption stops when the saturation value is 
reached. \B\ and IFI denote the 3-dimensional volume of I SI and the 2-dimensional surface area 



of F, respectively. In (2.11 ) W^ has to be understood as the trace of the cytosolic GDP-GTPase 



mol 



concentration and has units ^^ . The units of the various coefficients are given by 

m^ o J 

[D] = [du] = K] = — ; [66] = ^^; [6_6] = -• 
s mol • s s 

2.2.5. Reaction, Diffusion, and attachment/ detachment. Taking reaction and diffusion into ac- 
count, one obtains the following model, which is the basis of further considerations 

dtV = DAV in 5x1, (2.12) 

dtu = duATU + kivgQ[l- \ ] + k2V ^ ° h^ — — - on F x I, (2.13) 

\ I + A5U/ 1 + K^u u + k^ 

dtv = d,Arv-kivgo(l--^^)-k2V-^^^ + ks^^ + q on F x I (2.14) 
\ 1 + K^uJ 1 + K^u u + ki 
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with the flux conditions (2.10), (2.11). The model satisfies conservation of GTPase in the form 

— [ / Vdx+ {u + v)da{xU =0, 

where da{x) denotes integration with respect to the surface area measure. This equation confirms 
that the total number of cytosolic inactive plus membrane-bound active and inactive GTPase is 
constant over time. 

2.2.6. Non-Dimensionalization. For x £ B and t G I, we introduce non-dimensional coordinates 

where R > denotes a typical length, e.g. half the diameter of the cell. We represent this 
length as -R = ^/^I with 7 > and I = Im denoting the unit length. Furthermore, we use a 
dimensionless time 

This leads to transformed domain B := {^ £ W^ : RS^ £ B}, T := dB and time interval I := 
[0, -niT]. Non-dimensional Rab concentrations are defined through 

V := V, u := u, v := v. 

Cmax Cjnax Cjnax 

We introduce dimensionless quantities 

I\ 1 k2 I2 j.^ 
ai := —kigo, 02 := — , as := — ai, 04 := fes, 05 := , 

IX \B\ I%.e . d, ~ D 

06 := -^Cmaxppq^, a-Q := — — , d:=—, V := —. 

Note that tU^ is scale invariant in the sense that multiplying the system size by a constant factor 
a > does not affect this value. In particular all above constants are independent of the system 
size, which is solely represented by the dimensionless quantity 7. With these definitions, one 
easily verifies 

drV = DA^V in B x I, (2.15) 

drU = ApU + j{ I oi -|- (03 — oi) ) V — 04 ) on F x /, (2.16) 

V V a2 + uj a^ + uJ 

drV = dAf.v -I- 7( — ( ai -I- (03 — ai) ) v + 04 h g ) on f x / (2.17) 

V \ a2 + uj a^ + u J 

with the flux condition 

-DV^V -D = 'yq on f x / 

with 

q = aQV{l- {u + v))+-a.QV. (2.18) 

2.2.7. Reduction. We further reduce the non-dimensional model of the previous section. Our 
reduction is motivated by the observation that the cytosolic diffusion coefficient is much larger 
than that of the lateral diffusion on the membrane [28]. We thus assume V to be spatially 
constant, i.e. V = V{t) depends only on time but not on the ^ variable. If we also assume that 
the initial concentration of cytosolic GTPase is spatially homogeneous the concentration is then 
for positive times determined by GTPase conservation, 

V{T) = Vo-cj{u + v){^,T)da{0 (2.19) 
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where c := \B\ and where Vq is given by the initial conditions, 



In particular Vq = V'(O) if initially no membrane-bound GTPase was present. In the following, 



flux (2.18) and the conservation law (2.19). 



we then consider the system (2.16)-(2.17) of reaction diffusion equations on F x / including the 



The fully coupled system converges in the limit D — )• oo to this reduced model. However, no 
estimates for the difference between solutions to the respective models are at present available. 
The ratio between cytosolic and lateral membrane diffusion reported in the literature [28j is of 
order 10^. Nuinerical experiments for the full system with diffusion coefficients d = D = 10^ 
showed qualitative agreement with the reduction. 

3. Turing pattern 

In this section we investigate the stability properties and the possibility of Turing-type pattern 
formation for the dimensionless reduced model derived above. In the following we drop all tildes 
and denote the space and time variables by x and t respectively. This yields the system 

dtu = Aru + -ff{u,v), (3.1) 

dtv = dArv + 7 (-/(n, v) + q{u + v, v, V[u + v])) (3.2) 

where 

u \ u 



f{u,v) 



ai + (as - ai) 



v — 04- 



a2 + uj a^ + u 

q{u + V, V, V) = 06^(1 - (n + w))+ - Q-qv, 

and where y[u + v] is the non-local functional 

V[u + v] = Vq — c [u + v) dcr(x), 

with Vo given. For convenience we also define 

g{u,v) = -f{u,v) + q{u + v,v,V[u + v]). 



(3.3) 
(3.4) 

(3.5) 



System (3.1), (3.2) has to be solved on F x /. Our particular interest is to understand the 
effect of the non-local term in system (3.1), (3.2) on the stability properties. We can not use a 



predefined set of 'realistic' parameter values: first our model is very general and applies to several 
specific cases with different set of parameters; second kinetic rates etc. are difficult to obtain 
experimentally. We therefore rather investigate whether in principle, i.e. for some parameter 
values, our model allows for stationary states, whether stationary states of substrate-depletion 
type exists, and whether Turing type instabilities are possible. It is quite difficult to guess 
parameters that allow for example for Turing pattern formation. We therefore derive conditions 
for the parameters that are sufficient to ensure certain behavior, in particular showing that the 
Turing space is not empty. With such a set of parameters identified it is possible to explore the 
boundaries of the Turing space and then compare whether the parameter ranges are reasonably 
close to available estimates for 'realistic' values. 

To start with our stability analysis we first observe that spatially homogeneous solutions of 



( |3.lD , ( |3.2[ ) satisfy the ODE system 

dtu 
dtv 

where 



l{-f{u,v) + qo{u + v,v)) 



(3.6) 
(3.7) 



qo{u + v,v) = a6(Vb - c\T\{u + v)) (1 - {u + v))^ - U-ev. 



(3.8) 
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We set go{u, v) = —f{u, v) + qo{u + v,v). The set of values for u, v described by 



A := {u, V > : u + V < min{l, m}}, 



m :- 



ciri 



(3.9) 



is an invariant region for (|3.6l), (3.7), i.e. if the initial data are in this set the solution does not 



leave it. In fact we observe that at the boundaries of A we obtain 

/(0,w) > 0, goiu.Q) > for alio < u,v < min{l,m}, 
f{u, v) + gQ{u., v) < Q for all u, u > 0, n + i; = min{l, tti}. 

By these inequalities the conclusion follows. 

Under suitable conditions on the data we next show the existence of a stationary spatially 
homogeneous state («*,!'*) G A for (3.1), (3.2). This means that (tt*,t'*) has to satisfy 

= /(«*,«*), 
= go{u^,v^). 

The first equation is satisfied if and only if 

aiu{a2 + u) 



V = v[u] :- 



(as + u){aia2 + asu) 



We compute 



If we assume 



'j'[u]=0 <^ (aia2 + 03(05 — a2))ti +2010205^ + 010205 = 0. 



02 > 05, 
2oi02 < 03(02 — 05) 

we find that v[-] has a unique positive stationary point uq, 

_a/(o3 - ai)(a2 - 05) 



Uq 



In particular we have 



01O2O5 



03(02 - 05) - 0102 



+ 02^/0105- 



03(02 -05) -0102 



v'[u] < for all u > uq. 



By (3.14), (3.15) we estimate 

f 2,/E^ 1 s^ 

uo < Vaia5a2 — 7 r + 2 — 

^03(02-05) V03(02-05)>' 

V«l0502 



< 4- 



y/03(02 -05) 

For the maximum value vq := ^["00] we obtain 



Vo = 04 



02 + 2no 



01O2 + 03(05 + 2uo) 



Using ( 3.14 ) and ( 3.15 ) a short calculation yields 

02 04 



^'0 



< 



Vo > 



0102 + 0305 
0405 

0102 + 0305 



< 



> 



0204 
0305' 

0405 

0302' 



(3.10) 
(3.11) 



(3.12) 

(3.13) 

(3.14) 
(3.15) 



(3.16) 



(3.17) 



(3.18) 
(3.19) 
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If we then choose 

< -mmjm, 1|, 

03^5 4 

0305 

we deduce by (3.17), (3.18) that uq + vq < ^m.m{m,l}. 

In order to satisfy (3.11 ) we need to find u > with u + v[u] < min{rn-, 1} such that 



(3.20) 
(3.21) 



We evaluate 



and assuming 



«6 (^0 — c|r|(n + I'M)) (1 — (n + v[u])) — a-Qv[u] =: ^(u). 



^{uo) > —Vo-a-e 

4 0302 



Q-6 ^ Vo a3«2 



(3.22) 
ag 4 0405 

we see that ^{uq) > 0. On the other hand there exists ui > uq such that ui + v[ui] = min{?n,, 1} 
and we observe that 

$(ui) < 0. 

Since ^ is continuous we obtain from the intermediate- value Theorem that there exists uq < u^, < 
ui such that $(ii*) = 0. But this implies that {u^,v^) £ A, v^: = v[u^], is a stationary point of 
(3. 6), (3. 7). In summary we have proved the following Proposition. 



(3.23) 
(3.24) 
(3.25) 



p 


roposition 


3.1. Assume that the conditions 










a2 


> 


as, 










40204 


< 


asasminjm, 1}, 










4a4a5«-6 


< 


^0020306, 










ai 


< 


. r 03(02 -05) 
mm<^ — ^ 

I 2a2 


2" 


_ga|(a2 -05)^ 



mm 



{m,l})2} 



(3.26) 



are satisfied. Then there exists a stationary spatially homogeneous solution {u^,,v^,) E A of (3.1), 
(Q. 

The stationary point (it*,f*) is under suitable assumptions on the data linearly stable against 
spatially homogeneous perturbations. For a brief summary of the classical stability analysis and 
of the Turing mechanism for two- variable reaction-diffusion systems we refer to Appendix lAl 



Proposition 3.2. Assume that (3.23)-(3.26) hold and that moreover 

204(02 - 05) < asal, 

a_6 < a6c|r||l 



m\ 



(3.27) 
(3.28) 



are satisfied. Then (ti*,^*) is a stable stationary point of (3.6), (3.7). This system is in {u^,,v^,) 
of activator-substrate-depletion type, where u acts as an activator and v as substrate. 



Proof. We show that the stability conditions (A.5), (A. 6) are satisfied. We first observe that 
since oi < ^ by ( |3.26D 

u 



dvf{u,v) 



= ai + (03 
> 



01 , — 

for all u > 0. 



(3.29) 
(3.30) 
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For the function v[-] defined in (3.12) we have f{u,v[u]) = 0. This yields 
= duf{u,v[u]) = {duf){u,v[u]) + {d^f){u,v[u])v'[u] 



Since u* > no we deduce from (3.16) that w'[u*] < and obtain 

duf{u^,v^) = -dvf{u*.,v^)v'[u^] > 0. 

Furthermore we have 

duqo{u,v) = -aQc\T\[l + m - 2{u + v)) < 0, 
dyqo{u,v) = — a6c|r|(l + m — 2(u + w)) — a_6 < 0. 



(3.31) 

(3.32) 
(3.33) 



By (3.30)-(3.33) the stationary point (n*,f*) is of activator-substrate-depletion type. To check 
the criteria for Turing type instabilities we need to estimate combinations of derivatives. We first 
compute 



r. r 02(03 -ai) 

duf = — ZT^V 



0405 



(a2 + M)2 " (a5 + u)2' 

Evaluating this expression at {u,v) = {u^,,v^,) and using ( 3.12[ ) we deduce 

02(03 — ai)a4U 0405 

(02 + u)(a5 + 'u)(oia2 + 03U) (as + u)^ 
02O304 0405 



duf 



< 



(02 + u)(a5 +u)a3 (a5+u)2 
04(02 — 05)^ 



(02 + u){a5 + ti)2' 
We thus obtain in {u,v) = (n*,f*) that 

04(02 — 05)u 



(3.34) 



duf + d^go < 



oi + (03 -oi) 



u 



< 



a2 + n)(o5 + u)2 V ^ ' " " 02 + u 
— 06c|r| (1 + m — 2('u + v)) — o_6 
04(02 — 05)^ u 



(02 + n)(o5 + u)^ 



03 



02 + ti 



< 



(3.35) 



by (3.27). This verifies condition (A.5). We next estimate in {u,v) 

dufdvgo - dyfdugo 
= duf{-dvf + d^qo) - dyf{-duf + duqo) 
= duqoiduf - dvf) - a^eduf 

^ mi /I , -,/ , Wf 04(02 -05)^ 

> — oecirj [I + m — 2[u + v)) ' 
u 



{u^-,v^} 



u 



03- 



0-6- 



04(02 — osjn 



(02 +n)(o5 + n)^ 



02 +u)(05 +-0)2 a2+uJ (02 + u)(05 +m)2 

06c|r| (1 + m - 2{u + v)) (04(02 - 05) - 03(05 + u)^) + o_6a4(o2 - 05 

(3.36) 



By (3.27) and since l + 'm — 2{u-\-v) > |1 — r?T,| the term in the brackets in the last line is estimated 
by 

06c|r| (1 + m - 2(n + v)) (04(02 - 05) - 03(05 + n)^) + 0-604(02 - 05) 
< — 06c|r||l — m|o4(o2 — 05) + 0-604(02 — 05) < 0, 



where the last estimate follows from (3.28). Together with (3.36) the last equation implies 

dufdvgo - dyfdugo > 0, 



and therefore (A.6) holds. Thus the ODE system is linearly stable. 



D 
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We next evaluate the response of the fuh reaction-diffusion system to perturbations of the spa- 
tiaUy homogeneous solution (u^,, v^,) in direction of arbitrary smooth functions 99, V' : T x (0, T) — )• 
M. In particular we have to linearize the non-local functional V = V[u + v] that occurs in the 



source term q in (3.2). With this aim we consider a variation {us,Vs) of {u^,,v^,) in direction of 



Us,Vs : Tx (0,r) 



SG(-1,1) 



ls=0 



^*' '^sU=o 



V*, 



d_ 
ds 



Uo 



"P, 



d_ 
ds 



s=0 



i,. 



The corresponding linearization of V is then given by 



d I 



V{us 



A r 



(V7 + ^)dc7(e). 



(3.37) 



For the linearization of (3.1), (3.2) we therefore obtain 



(3.38) 
(3.39) 



dtip = dAri) + 'y[- duf{u^,v^)<p - dvf{u^,v^)i' + diq{u^ + v^:,v^:,V^){(p + 11;)) + 
+ 'y(d2q{u^ +v^,v^,V^)ip - cd^q{u^ +f*,f*,T4) / (c^^ + V") dcr(^) j. 

We next decompose the direction of perturbation (c^, ^) in L'^iT) in a part that is spatially homo- 
geneous and a part that is orthogonal to the constants. Since spatially homogeneous perturbations 



have already been analyzed in Proposition 3.2 it suffices to assume that we have 



ipdLa{x) 



iljda{x) 



0. 



Equation (3.37) shows that for such directions V is unchanged to first order. Therefore, with re- 



spect to variations in direction of functions that are orthogonal to the constants, the linearizations 



of (3.1), (3.2) coincide with that of the system 



dtu = Aru + ^f{u,v), 

dtv = dArv + 7 {-f{u, v) + qi{u + v, v)) 

in (u*, v^k), where 

qi{u + v,v) = q{u + v, v, K) = a6K(l - (u + v)) - a-ev, 

K = V[u^: + w*] For convenience we define 

gi{u,v) = -f{u,v) + qi{u + v,v). 



(3.40) 
(3.41) 



(3.42) 



Thus we see that the non-local term in the full system ({S.lh, (3.2) leads to a difference in the 



linearization with respect to homogeneous or heterogeneous perturbations, that can be understood 
as a change (from qo to qi) in the source term. Below we show that a range of parameters exists 
where {u^,v^) is an unstable stationary state. We first need some estimates on m^,,^^, to prepare 
the proof. 



Lemma 3.3. Assume (3.23), (3.24), (3.28), and 



aia2 < 



«3 



1 + 4- 



(3.43) 
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0402 



0305 
u^: > -min{l,m}, 



w* > 



04(02 + l)(fl3 + 1) min{l,m} 
2(05 + I)(a2 + 2)03 



(3.44) 
(3.45) 
(3.46) 



holds. 



Proof. By (3.18) and v^: < vq we deduce (3.44). From 
solution of 



{u^,v^:) = we obtain that {u^,v^) is a 







06c|r| [m — {u + v)) (1 — (u + f )) — a-QV. 



(3.47) 



If we denote by u[v] G (0,niin{l,r?T,}) the solution of (3.47) for given < v < minjl,???-} we 
observe that u is decreasing in v. By (3. 24), (3. 44) we have v^: < vi := ^ minjl, m} and therefore 



u* > ui, ui := u[vi]. By (3.47) and (3.28) we get 



Ul 



hm+1- 2t;i) - ]{m + l- 2viY -{m- vi){l - vi) + ^^ 
2 V 4 flee r 



> -(2niax{l, m} + min{l,m}) — \ -rifn ~1)^ + t|1~'^| niin{l,m} 
min{l,m}, 



which proves (3.45). Next we obtain from f(u,v) = for {u,v) = (ti*,f*) that 



04(02 + l)(a3 + l)n 
(05 + n)(aia2 + 03^*) ~ (05 + l)(aia2 + 03^) ~ (05 + l)(a3 + 2)03 ' 



aiu{a2 + u) 



> 



0.4(0.2 + l)^t 



> 



where we have used (3.23) and (3.43). By (3.45) this yields ( 3.46| ) 



D 



Theorem 3.4. Let («=„, v^,) be the stationary state found in Proposition 3.1 and let the parameters 



in (3.1), ( |3.2[ ) satisfy all the conditions (3.23)-(3.28). If in addition 
■min{l,7TT,}a3 min{l,m}^a3(a2 — 05) 



ai < mm 



.{: 



1) /' 



d > 



d > 



2o2(a2 + l) ' 4a2(a2 + l)2(a2 + 

2(03(0! + 2) + (a2 + l)(ai + l){aGVo + a-e)){ol + 2)(a5 + 1)^ 

min{l, m}a3a4(a2 — a5)(a3 + 1) 
A{al + 2)2(a5 + lf{alai{a2 - 05) min{l, m} + 4(a| + 2)05^0(02 + l)(a5 + 1)^) 
03(03 + 1)04(02 — 05)2 min{l, m}2 



(3.48) 
(3.49) 
(3.50) 



then there exists 7 > such that {u^:,v^:) is linearly unstable. 



Proof. We show that the conditions (A. 13), (A. 14) stated in Appendix [A| for the instability of 
(3.40), (3.41 ) are satisfied. Let us start with (A. 13) by estimating the different partial derivatives. 



duf{u*,v^) 



02(03 — 01)0414(05 + u) — 0405(0102 + 03U)(02 + u) 
(02 + ti)(05 + «)2(oia2 + 03«) 

04(a3(«2 - a^)u'^ - 0102(205^ + u^ + 0205)) 

(02 + u)(05 + ti)2(oi02 + 03«) 
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We next observe that by ( |3.23D , ( |3.48D and ( |3.45D 

0102(205^ + 14^ + 0205) < 0102(1 + 02)^ < 

as 



13 



as 



1 + 0: 



2(02 -05JU 



0102 < 



1 + oi 



(3.51) 



We therefore deduce that 



duf{u*,v*) > 



04(03 



03 



l+a 



7)(o2 - a^)v? 



(02 + n)(a5 + m)2(^ + 03)14 



040|(02 — 05)ti 



> 



(02 + l)(05 + l)2(02+2) 

0403(02 — 05) minjl, m} 

2(02 + l)(05 + l)2(0| + 2)' 



(3.52) 



Next we estimate 



ai + (03 - oi) 



> 



> 



0102 + 03U 
02 + u 
as(ai + 2) 

'(02 + l)(0i + l) 



02 + U 

(oeVb + a_6) 



Og 



Vo-c|r|(u + v))+o_6) 



(o6T4) + a_6), 



where we have used (3.51). Together with (3.52) we deduce from (3.49) that (A. 13) holds. 
Next we verify the condition (A. 14), i.e. 

D := {dduf + d^gif - Ad{dufd,gi - d^fduQi) > 0. 

We compute 



dufdvQi - dyfduQi = dufdvQi - dvfduQi 
and obtain for the left hand side 

D = (f{duff + 2d( - dufidj + d^qi) + 2dyfduqi) + (d^gif 
> d 



Moreover 



dvf 



d{dufY-2[dufd,f-2d,fduqi 

03(0^ + 2)u 



O1O2 + 0314 



< 



02 + 14 (o2+n)(a| + l) 

oeK < oeVb. 



< 



a3(a3 + 2) 
(o2 + l)(oi + l)' 



This yields 



D > d{duff 



d 



2o3(o| + 2) 



4o3(o| + 2)o6Vb 



(02 + l)(oi + l){duf) (02 + l)(oi + l){dufY- 

We then deduce (A. 14) from ( 3.50| ) and (3.52). The conclusion now follows from (15]. D 

The above conditions ensure that perturbations with eigenvectors of the Laplace-Beltrami oper- 
ator on r corresponding to eigenvalues in a certain interval are unstable. As this interval scales 
linearly with 7 we obtain a range of values for 7 where a Turing instability exists. 

We conclude this section with some comments on the implications of the conditions derived 
above. 



14 



A. RATZ AND M. ROGER 



Remark 3.5. By Theorem 3.4 parameters satisfying ( 3.23 )-( 3.28) and ( 3.48 )-( 3.50) belong to 



the Turing space where diffusive instabihties exist. Some of these conditions can be easily in- 
terpreted. The requirement 02 > 05 concerns the Michaelis-Menten constants appearing in the 
catalyzed reactions: the affinity of GEF towards activated GTPase (forming the GEF-GTP- 
GTPase-effector complex) has to be higher than the affinity of GAP towards activated GTPase. 
Several conditions require oi to be (much) smaller than 03, which means that activation by the 
GEF-GTP-GTPase-effector complex is stronger than activation by single GEF molecules, which 
in fact has been reported for example in the case of Rab5 GTPase [20]. The conditions (3.49), 



(3.50) for a Turing instability are most critical, as a substantially larger lateral diffusion coefficient 



for inactive GTPase compared to the lateral diffusion coefficient for active GTPase is required. 
We investigate below whether this condition is due to the particular choices of kinetic and sorption 
laws or rather a general feature of the kind of (reduced) model that we are considering. Finally, 
the condition on 7 requires a certain minimal size of the system to allow for a Turing instability. 

4. Stability analysis for equal lateral diffusion 

As in most applications no substantial difference in the diffusion coefficients of the GDP-bound 
and GTP-bound GTPase is present we investigate in this section the possibility of Turing pattern 

1 of equal lateral diffusion. The non-locality of our model - the 

changes the classical stability analysis. 
1 might become possible. 



in (3.1), (3.2) for the case d 



remnant of the full 2D-3D coupling in our reduction 

Therefore, in contrast to the classical case, Turing pattern for d 

However, we show here that in our simple model this is not the case. 

The set-up in this section is as follows. We assume a system of the general form 

dtu = Aru + jf{u,v), (4.1) 

dtv = Arw + 7 (-/(u, v) + q{u + v, v, V[u + v])) (4.2) 

where u, v denote GTP-bound and GDP-bound GTPase concentrations, respectively, and where 
y[M + w] represents the cytosolic (inactive) GTPase concentration that is given by the mass 



conservation condition (3.5). The nonlinear terms f,q account for the activation/deactivation 



processes and from the attachment/detachment of GTPase at the membrane. For q we assume 
that 

diq < 0, d2q < 0, dsq > 0, (4.3) 

which are natural condition with respect to the interpretation of q as the flux induced by ad- and 



desorption of GTPase at the membrane. As before, the system (3.1), (3.2) has to be solved on 
Fx /. 



We assume a spatially homogeneous stationary point {u^,v^) of the ODE reduction of (4.1), 



(4.2) 



dtu = -ff{u,v), 

dtv = ^ {-f{u,v) + qo{u + v,v)) 
where 

qo{u + v,v) = q{u + v,v,Vo{u + v)), Vo{u + v) = 
We set as before 

g{u,v) := -f{u,v)+q{u + v,v,V[u + v]), go{u,v) = 
The main result of this section is the following theorem. 

Theorem 4.1. Assume that (ti*,^*) is of activator-substrate depletion type, that is 
duf{u^,v*) >0 dyf{u^,v^) > 0, 

dugo{u^,v^) <0 dygo{u^,v^) > 0. 



(4.4) 
(4.5) 



Vo-c\T\{u + v). 
-f{u,v) + qQ{u + v,v) 



(4.6) 
(4.7) 
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Then no Turing type instability of (4.1), (4.2) exists. 



Proof. The conditions (A. 5), (A. 6) to ensure that (u*,f=i,) is a stable stationary point of (4.4), 
(4.5) are 

> duf + d,go = duf - d,f + diq + d2q + d^qV^, (4.8) 

< a„/ • d^go - d,f ■ dugo = duf{diq + d2q + dsqV^) - dj{diq + d^V:,). (4.9) 



This corresponds to the stabihty of (4.1), (4.2) in (u*,^*) with respect to spatially homogeneous 
perturbations. Therefore, for the instabihty with respect to general perturbations it is sufficient 
to consider perturbations in direction of functions perpendicular to the constants. As above we 
observe that such perturbations leave V[u^ v] unchanged. The respective linearization corresponds 
to that of the system 

dtu = ^TU + -ff{u,v), (4.10) 

dtv = dArv + 7 {-fiu, v) + qi{u + v, v)) (4-11) 

in (u*, V*), where 

qi{u + v,v) = q{u + v,v,V^.), V^ =V[u^ + v^,] = Vo{u^,v^). (4.12) 

We set as before gi{u,v) = —f{u,v) + qi{u + v,v). Then the conditions (A. 13), (A. 14) for the 
instability of (4.10), (4.11) with respect to spatially heterogeneous perturbations yield 

< duf + d,gi = duf - dj + diq + d2q, (4.13) 

< {duf + d^gif - 4{duf ■ d^gi - dyf ■ dugi) 

= {duf - d,f - diq + d2qf - ^duf ■ d2q + 4diq • d2q. (4.14) 



By our assumptions (4.6), (4.3) we observe that the last condition is automatically satisfied. On 



the other hand, we obtain from ( |4.9| ) that 

< {duf - d,f){diq + d^v^) + dufd2q. 



(4.15) 



By (4.13) and (4.3) we deduce that duf — d^f = —{diq + d2q) > 0- Using again (4.3) and using 
that Vq <Q this yields that the first term on the right-hand side in (4.15) is nonpositive. But we 



also find by (4.6), (4.3) that the second term is nonpositive. This is a contradiction. Therefore 
the system has no Turing-type instabilities. D 



5. Numerical Approach 

We use a finite element discretization similar to the one described in [19]. It is implemented 
in the adaptive finite element toolbox AMDiS [32j. 

5.1. Discretization. Following the surface finite element method described in [6j, we choose a 
triangulated discrete approximation T/j of the membrane V and a triangulation Th- We split the 
time interval [0,T] by discrete time instants to < ti < • • • < t^/, from which one gets the time 
steps Aim := tm+i — ^m, rn, = 0, 1, . . . , M — 1. Given initial conditions «(-, 0) = uq, f (•, 0) = vq 
with uq, vq S H^{T}i) and time discrete solutions u^"^', v^"^' G H^{Tii), m = 1, . . . , M, we linearize 
all nonlinear terms 

/•(^(m+l) -^("^+1)) ?« f(ti(™) 7)(™)l ■ " '" '""-^ ''"'"'' '^ ~ "" 



,v-->) + Vf{u^'"^>,v^"'>) 



(m) „,(m)\ 



yim+1) _ ^(m) 



and 






+V(„,,)9(n(-),^(-),y(-))- 



.(m+l) 



(m) 



16 A. RATZ AND M. ROGER 

We introduce test functions r]^,7]'" G H^(rh) and end up with a weak formulation and semi- 
implicit time discretization for u^"i+i) ^ ^{m+i) ^ H^(rh) of ([sl]), ^^ 



At 






// (rn+l)\ 



— / t.(™)r/^- / Fe(n(™),f(™))?7^+ / Qe(w^'"\ v^"\ "^^™^)??'' Vr/'' G //^(r,,), 



Aim JTh JFh ^Th 

where 

((m) 

Furthermore, the non-local relation for V^"^> is treated explicitly by 

for given Vq and the inner Bh of F/i. To discretize in space, let Yh the finite element space of 
globally continuous, piecewise linear elements. In addition, with {ipi)i the standard nodal basis 

of V/i and u^™" ,Vf^ G V/i we write uj^ = Yl^i Vi and v^ = Yl^i vi with 

i i 

^ira+i) ^y{m+i) ^ ^_ Furthermore, we define C/('"+i) = {ut^\ and ^("^+1) = (1//™+'^)^. This 
leads to the linear system of equations 

1 

At 






At. 



l__/v/\/('") — i^cxpl _|_ QGxpl 



with 



A = {Aij ) Aij = ( VV'i , ViPj )r^ , 

A2 = (4) 4 = (^(V^l™))^., VV,)r„ 

FT"' = i{Fr%) {Fr"% = id.f{ui"'\vlr^)^,,^,)r,, 
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where (•, •)r^ denotes L^-scalar product. The above hnear system has to be solved in every time 
step, which is done by a stabihzed bi-conjugate gradient method (BiCGStab). 

5.2. Numerical Results. First, we present numerical results reproducing the results of the 
stability analysis in Section [3| To be more precise, we choose a set of paraineters fulfilling the 
conditions of Theorem |3.4| sufficient for instability, which is the basis of our further numerical 
investigations: 

d = 1000; ai = 0; 02 = 20; 03 = 160; 04 = 1; 05 = 0.5; a^ = 0.1; a_6 = 1; 7 = 400. (5.1) 

Furthermore, we consider the unit-sphere T = S"^ and random initial conditions uq,vq : Th — )• 
[0,0.02] and Vq = 10. Fig. [2] shows the corresponding discrete solutions Uh,Vh at different times. 
A stationary pattern with a single spot appears. 







V 

1,04 

;o.Q3 
;o,oi 

Q 



Figure 2. From left to right: the discrete solutions Uh, Vh for t 
t = 5, and t = 25. 



tQ = 0,t = 0.5, 



5.2.1. Varying Parameters. Based on the choice of parameters (5.1) we investigate the influence 



of varying parameters. First, we observe that doubling the parameter 02 leads to spatially homo- 
geneous stationary solutions, whereas halving 02 leads to stationary patterns with two maxima 
for Uh (see Fig. |3]). 







0,026148 
[0,024 



0,012 
0,011547 



Figure 3. From left to right: the discrete solutions Uh,Vh for 02 = 40 (left) and 
02 = 10 (right) t = 25. 



Additionally, we observe that halving the parameter 03 leads to spatially homogeneous sta- 
tionary solutions, whereas doubling 03 leads to stationary patterns with two maxima for u^ (see 
Fig. 11). 

Another interesting behavior can be seen by varying the diffusion constant d. While the esti- 



mates (3.49) and (3.50) lead to a condition d > 790 sufficient for instability, an exact computation 

; equality in one of relatio 
100 and instability for d 



yields a maximal diffusion constant dc ~ 101 satisfying equality in one of relations (|A.13|), (|A.14|). 
This is reproduced in Fig. [5] showing stability for d 



105. 
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1,022472 
0,02 
0,0175 
0,015 
0,0126 
0,01 
0,008214 



Figure 4. From left to right: the discrete solutions Uh,Vh for 03 = 80 (left) and 
as = 320 (right) t = 25. 






0,103969 
0,1 



0,08 
0,06 



-0,04 
0,035732 



Figure 5. From left to right: the discrete solutions Uh,Vh for d = 100 (left) and 
d = 105 (right) t = 25. 



5.3. Comparison with 'realistic' parameter ranges. A full set of realistic parameters is 
not available, but there are several in vivo and in vitro measurements giving some estimates or 
average values. For the case of the Cdc42 GTPase cycle in yeast cells we have collected data from 
[9], [13], [E], and [7] and evaluated our model for the following set of values, 

-5 



k. 



ki := 1.056831769 • 10~ , 

fc4 := 18.92448788, k^ := 0.1056831769 

b-e ■= 0.133, go := 37848.97575, d„ : 



0.1056831769-10" 



10" 
= 2.5-10 



ks := 946.2243938 
k^5 ■= 0.3, be := 0.3170495307 
max := 47311.21969, 



10" 



'15 



Vo := 4.894264108 - 10 
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where the units are as given in Section [2j 

We find for these values that a homogeneous stationary state of activator-substrate-depletion 
type in fact exists. For a Turing-type instability we need a lateral diffusion for activated GTPase 
of order 10~'^m^s"^, which is unrealistically large. Nevertheless assuming such a value we obtain 
a condition on the spatial scale: We find that R has to be at least of order 10~^m, which is close 
to the typical diameter of a yeast cell. We therefore see that the critical condition is in fact the 
large difference in diffusion. 

The uncertainty in the parameter values is quite large: there are not enough data for individual 
proteins available and the parameters chosen above are therefore a composite of available data. 
Many measurements are taken in vitro and are only an estimate for in vivo conditions. Comparing 
different sources and different GTPases one may find variations up to order 10 in the parameters. 
Within the set of parameter values allowing for Turing instabilities we find values that are within 
the range of realistic values, except that the ratio d of lateral diffusion values has to be of order 
at least 10^. Such a large difference in free lateral diffusion for active and inactive GTPase seems 
unrealistic. However, heterogeneities of the cell membrane may lead to differences in the 'effective' 
diffusion speed, see the discussion below. 

6. Discussion 

The GTPase cycle presents an example for a coupled system of processes in the inner volume 
and on the outer membrane. We have proposed a mathematical model in the form of a fully 
coupled PDE system. A two- variable reduction yields a non-local reaction-diffusion system on 
the membrane. With the interest in finding mechanisms that support the emergence of cell 
polarity we have investigated pattern forming properties. We have shown that the reduced model 
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in principle supports Turing type diffusive instabilities, but - as for general local RD systems 
- needs large differences in diffusion constants. In numerical simulations we have confirmed 
our theoretical findings and have explored the type of pattern produced and the infiuence of 
parameter changes. Both the formulation of the model and the numerical schemes are prepared 
to investigate more complex models and to incorporate additional features. 

In similar but different models diffusive instabilities have been shown to exist 12], j9], or not 
to exist [I] , [3] (unless a phase separation force is added) . Particularly interesting is a comparison 
with the work of Goryachev and Pokhilko |^ . Their model is more detailed in the set of variables 
they consider and also accounts - at least partially ~ for the different dimensionalities of the 
processes: the membrane is thought as a thin compartment with positive volume and the cell 
as an adjacent bigger compartment. Concentrations on the membrane and in the inner cell are 
weighed with a factor that accounts for the different sizes of the volumes. A major difference 
to our work is that the mathematical analysis in [9| treats all variables on one common domain 
of definition (for both cytosolic and membrane-bound quantities). In our approach on the other 
hand we distinguish explicitly between the cytosolic GTPase variable V defined in B and the 
membrane variables defined on T, which then makes laws for fluxes from the cytosol to the 
membrane necessary and meaningful. Nevertheless there are more similarities between these two 
approaches. Also in [9] a non-local reduction to a two-variable system is given. The substrate 
there is represented by the sum of cytosolic and membrane bound GTPase and inherits a higher 
diffusion constant than the solely membrane bound active form. In view of pattern forming 
properties their model then produces Turing instabilities, in more realistic parameter ranges than 
our model. 

One of the main features of our reduced model is that large differences in the lateral diffusion 
constant for active and inactive GTPase are required to obtain a Turing instability. Mathematical 
models with a simpler (but more artiflcial) dimensional coupling on the other hand do support 
diffusive instabilities. However, in these studies differences in lateral diffusion were put more 
directly into the model and Turing instabilities then appear as a consequence of this model design. 
Our analysis demonstrates that it is not clear whether large cytosolic diffusion is sufficient for 
cell polarization by a Turing mechanism. To clarify this more detailed studies are necessary. The 
absence of realistic Turing patterns in our model might be a consequence of our reduction to the 
infinite cytosolic diffusion limit. This leads to constant cytosolic GTPase concentration, whereas 
the Turing mechanism intimately relies on spatial heterogeneity. We have performed numerical 
experiments to compare the behavior of the fully coupled system and the reduced model. For 
large but finite cytosolic diffusion our simulations for the full model agreed qualitatively with 
corresponding simulations for the reduction. However, it is difficult to numerically explore the 
Turing space of the fully coupled model and error estimates for differences of solutions to the 
reduced model and the three-variable system are at present not available. 

In our model various properties of the GTPase cycle in living cells have been neglected that 
in principle may contribute to polarization. In a recent paper by Butler and Goldenfeld [1] 
it was shown that the inclusion of intrinsic noise in reaction-diffusion systems can lead to the 
formation of Turing-type 'quasipatterns' for parameter values substantially different from the 
Turing space for the corresponding noise-free system. The heterogeneity of the plasma membrane 
might influence polarization in living cells: Microdomains with different lipid composition are 
present, and active and inactive forms of GTPase possibly associate with different preferences to 
distinct microdomains. For the case of Ras GTPase in [21| a severe reduction in Ras mobility 
has been observed upon activation. This effect might introduce a difference in diffusion sufficient 
for Turing patterns. Finally, other mechanisms different from Turing pattern formation might be 
responsible for cell polarization. Possible candidates are a wave pinning mechanism j22j . or an 
intrinsically stochastic mechanism that relies on a particle based approach (as proposed in [1]) 
instead of a continuous model. 
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A better understanding of symmetry breaking properties in models where systems of different 
dimensionalities are coupled remains important. We have proposed a more detailed coupling 
model that presents a good basis for future extensions. A stability analysis of the fully coupled 
model and the inclusion of possible additional contributions to cell polarization are interesting 
tasks for further studies that might lead to a more complete picture of the origin of cell polariza- 
tion. 



Appendix A. Turing instability 

Since in our GTPase cycle model the classical conditions for Turing type pattern have to 
modified by the non-locality of the source term, we briefly outline the analysis of the classical 
Turing mechanism. We follow here |15' Section 5.3]. Let a spatial domain fl C M" be given and 
consider a system of reaction-diffusion equations 

dtu = Au + -ff{u,v), (A.l) 

dtv = dAv + -fg{u,v), (A. 2) 

where u = u{x,t), v = v{x,t), d > 1, 7 > and where / : M^ — t- M, </ : M^ — t- M are given 
functions. We complement ( A.l[ ), (A.2) first by initial conditions 

n(x,0) = n°(x), v{x,0) = ^^(x), 

for given u^,v^ : — )• M, and second by zero Neumann-boundary data 

Vu • vqu = Vv ■ I'n = on 90, 

where i^n denotes the outer normal of il.. With this boundary conditions the following analysis 
carries immediately over to the case that the spatial domain is given by a compact closed smooth 
hypersurface in M", with the only difference that the Laplace operator has to be replaced by the 
surface Laplace-Beltrami operator. 

The Turing mechanism is described by a stationary point (u* , v* ) that is spatially homogeneous 
and linearly stable under spatially homogeneous perturbation, but that is linearly unstable under 
heterogeneous perturbations. We therefore consider now (u*,f*) € M? with 

fiu^,v^) = 0, g{u^,v^) = 0. 

For spatially homogeneous solutions ( A.l[ ), (A.2) reduce to the ODE system 

dtu = lf{u,v), (A. 3) 

dtv = ig{u,v). (A. 4) 

The condition of linear stability then reduces to the conditions that the trace of the Jacobian 
D{f,g) is negative and that the determinant of D{f,g) is positive, i.e. 

duf{u*,v^) + d.ag{u^,v^) < 0, (A. 5) 



< 0, 

duf{u^,v^)d^g{u^,v^) - dyf{u^,v^)dug{u^,v^) > 0. 



(A.6) 



functions c^, 7/^ : ri x (0, T) 



The linearization of (A.l), (A.2) in («*,«*) for perturbations in direction of arbitrary smooth 



is given by the system 

+ 7 



1 
d 






(A.7) 



We next take the complete orthonormal basis {wj)j(z^^ of L'^{T) given by eigenvectors of the 
Laplacian with respect to zero Neumann boundary data. 



-Aw 



Vwi 



3 



XjWj 



in Q, 
on dQ, 



(A.8) 
(A.9) 
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for j = 0, 1, 2, . . . and where Aq = < Ai < A2 < • • • denote the corresponding eigenvalues. By 



the orthonormaUty condition and (A.8), (A. 9) we have 
ll'^illL2(n) = 1 



WjWi dx 



for aU i G No, 

/ ViUj • Vwi dx = for ah z, j G Nq, « / j. 

o Jn 

We then decompose ip{-,t),ip{-,t) with respect to this orthonornial basis, 

ip{x,t) = ao{t)wo + '^/3j{t)wj{x), 

ij{x,t) = Po{t)wo + ^Pjit)wj{x). 



(A.IO) 
(A.ll) 



Inserting this representation in (3.38), (3.39) and taking the L'^{T) scalar product with Wi by the 
orthonormality and (A.IO) the equation (A.7) decomposes into linear systems 

'1 0^ 



dt 



a 



Q, 



Pi) \0 d) \/3iJ \dug{u*,v^,) dyg{u^,v^)J \/3, 



duf{u*,v^) dvf{u^,v^) 



(A.12) 



for i = 0, 1, . . .. The case i = corresponds to the case of spatially homogeneous perturbations. 



In order to have an instability of (A.l), (A. 2) we therefore need that for an i G N (A.12) is 



unstable. This gives the necessary conditions [15, Theorem 5.3.1] 

dduf{u^,v^) + d,jg{u^,v^) > 0, 

(A.13) 

[dduf{u^,v^) +dvg{u*,v*)) - Ad(duf{u*,v*)dvg{u*,v^) - dvf{u*,v*)dug{u*,v*)) > 0. 

(A.14) 

In order to have an instability under these conditions it is sufficient that there exists an eigenvalue 
Aj, i G N, such that 

/x_ < Aj < /i+ 

where n = IJ,± are the roots of the quadratic equation 

df/ -7((i9„/(u*,f*) +5t,g(u*,w*))/i + 7^(9„/(n*,u*)9t,5'(n*,f*) - dvf{u^,v^)dug{u^,v^)) = 0. 
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